Part 1: Description

Step1: Downlaod the library

Note: install.packages(c(“tidycensus”, “tidyverse”, “sf”, “mapview”, “geepack”))

Step2: Get an API Key

Note: You only need to do this once per device using install = TRUE.

Step 3: Prepare census data

Step 4: Interactive map

Year: 2025

Year: 2020 to 2024

Year: 2020 to 2025

Part 2: Statistical analysis

Cross sectional analysis in 2025

## 
## Call:
## glm(formula = N_SHOOTINGS_2025 ~ N_BROKEN_LIGHTS_2025, family = poisson, 
##     data = final_cross_section_2025)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.91705    0.19716  -9.723  < 2e-16 ***
## N_BROKEN_LIGHTS_2025  0.04535    0.01503   3.018  0.00254 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 252.75  on 309  degrees of freedom
## Residual deviance: 244.93  on 308  degrees of freedom
## AIC: 365.33
## 
## Number of Fisher Scoring iterations: 6
##          (Intercept) N_BROKEN_LIGHTS_2025 
##            0.1470408            1.0463913

Longitudinal analysis during 2020 to 2025

## 
## Call:
## geepack::geeglm(formula = N_SHOOTINGS ~ N_BROKEN_LIGHTS + created_year, 
##     family = poisson, data = final_longitudinal_data, id = GEOID, 
##     corstr = "unstructured")
## 
##  Coefficients:
##                   Estimate    Std.err   Wald Pr(>|W|)    
## (Intercept)     448.829995  56.799488 62.442 2.78e-15 ***
## N_BROKEN_LIGHTS   0.012353   0.004656  7.039  0.00798 ** 
## created_year     -0.222397   0.028100 62.639 2.44e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = unstructured 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)    2.459  0.2596
##   Link = identity 
## 
## Estimated Correlation Parameters:
##           Estimate Std.err
## alpha.1:2   0.7595 0.09163
## alpha.1:3   0.4816 0.08234
## alpha.1:4   0.4816 0.07975
## alpha.1:5   0.5573 0.08891
## alpha.1:6   0.2690 0.06500
## alpha.2:3   0.6232 0.09262
## alpha.2:4   0.6866 0.08908
## alpha.2:5   0.6844 0.10104
## alpha.2:6   0.3365 0.08126
## alpha.3:4   0.4244 0.06458
## alpha.3:5   0.5023 0.09885
## alpha.3:6   0.2154 0.05049
## alpha.4:5   0.4427 0.05812
## alpha.4:6   0.2284 0.06247
## alpha.5:6   0.2165 0.05758
## Number of clusters:   310  Maximum cluster size: 6
##     (Intercept) N_BROKEN_LIGHTS    created_year 
##      8.402e+194       1.012e+00       8.006e-01

Part 3: Stratified analysis

Cross sectional analysis in 2025

Step 1: Prepare dataset

Step 2: see if there is association between broken street lamp and shooting case

Day time
final_cross_section_2025_DAY <- manhattan_pop |> 
  left_join(shootings_per_tract_25_DAY, by = "GEOID") |> 
  left_join(lights_per_tract_25, by = "GEOID") |> # lights_per_tract_25 calculated in Part 2
  mutate(
    N_BROKEN_LIGHTS_2025 = replace_na(N_BROKEN_LIGHTS_2025, 0),
    N_SHOOTINGS_2025 = replace_na(N_SHOOTINGS_2025, 0)
  ) |>
  st_drop_geometry()

poisson_model_2025_DAY <- glm(
  N_SHOOTINGS_2025 ~ N_BROKEN_LIGHTS_2025, 
  family = poisson, 
  data = final_cross_section_2025_DAY
)
summary(poisson_model_2025_DAY)
print(exp(coef(poisson_model_2025_DAY)))
Night time
final_cross_section_2025_NIGHT <- manhattan_pop |> 
  left_join(shootings_per_tract_25_NIGHT, by = "GEOID") |> 
  left_join(lights_per_tract_25, by = "GEOID") |> 
  mutate(
    N_BROKEN_LIGHTS_2025 = replace_na(N_BROKEN_LIGHTS_2025, 0),
    N_SHOOTINGS_2025 = replace_na(N_SHOOTINGS_2025, 0)
  ) |>
  st_drop_geometry()

poisson_model_2025_NIGHT <- glm(
  N_SHOOTINGS_2025 ~ N_BROKEN_LIGHTS_2025, 
  family = poisson, 
  data = final_cross_section_2025_NIGHT
)
summary(poisson_model_2025_NIGHT)
print(exp(coef(poisson_model_2025_NIGHT)))

Step 3: Compare through the table

Poisson Regression Results: Broken Street Lights vs. Shooting Cases (2025)
Time Period IRR (Incidence Rate Ratios) P-value
Daytime (6 AM - 7:59 PM) 1.057 0.0081
Nighttime (8 PM - 5:59 AM) 1.036 0.0973

Longitudinal sectional analysis during 2025

Step 1: Prepare the dataset

Step 2: See if there is association between broken street lamp and shooting case

Day time
poisson_gee_model_DAY <- geepack::geeglm(
  N_SHOOTINGS ~ N_BROKEN_LIGHTS + created_year, 
  data = final_longitudinal_data_DAY,
  id = GEOID, 
  family = poisson,
  corstr = "unstructured" 
)

summary(poisson_gee_model_DAY)
print(exp(coef(poisson_gee_model_DAY)))
Night time
poisson_gee_model_NIGHT <- geepack::geeglm(
  N_SHOOTINGS ~ N_BROKEN_LIGHTS + created_year, 
  data = final_longitudinal_data_NIGHT,
  id = GEOID, 
  family = poisson,
  corstr = "unstructured" 
)

summary(poisson_gee_model_NIGHT)
print(exp(coef(poisson_gee_model_NIGHT)))

Step 3: Compare through the table

GEE Poisson Regression Results: Broken Street Lights vs. Shooting Cases (2020-2025)
Time Period IRR (Incidence Rate Ratios) P-value (Wald Test)
Daytime (6 AM - 7:59 PM) 1.01 0.058
Nighttime (8 PM - 5:59 AM) 1.02 0.000